#ifndef __XMESH_DEC_H
#define __XMESH_DEC_H

#include "primitive.h"
/** 
 * @file dec.h
 * This file defines discrete exterior calculus on triangular mesh */

template <typename Value> class DecTriMeshVertex ; 
template <typename Value> class DecTriMeshFace ;

template <typename Value>
struct DecTriMeshPrimitive : public MeshPrimitive <Value>
{
//	typedef Value            ValueType ; 
	typedef DecTriMeshFace <Value>     FaceType ;
	typedef DecTriMeshVertex <Value>   VertexType  ;
//	typedef SubHalfEdge <Value> HalfEdgeType ; 
//	typedef Edge <Value>     EdgeType ;
	typedef Vector<Value>    VectorType ;
	typedef VectorType       PointType ;

};

template <typename Value>
class DecTriMeshVertex : public Vertex<Value>
{
protected :
	typedef typename DecTriMeshPrimitive<Value>::ValueType    ValueType; 
	typedef typename DecTriMeshPrimitive<Value>::FaceType     FaceType;
	typedef typename DecTriMeshPrimitive<Value>::VertexType   VertexType;
	typedef typename DecTriMeshPrimitive<Value>::HalfEdgeType HalfEdgeType ; 
	typedef typename DecTriMeshPrimitive<Value>::EdgeType     EdgeType ; 
	typedef typename DecTriMeshPrimitive<Value>::VectorType   VectorType ;
	typedef typename DecTriMeshPrimitive<Value>::PointType    PointType ;
	typedef typename Vertex<Value>::VVIter VVIter ;
public :
    /**
	 * Calculate the area of local averaging region. 
	 */
	const ValueType get_local_area ()  const 
	{
		ValueType area = ValueType(0);
		PointType center (ValueType(0), ValueType(0), ValueType(0));
		unsigned count = 0; 


		for (VVIter iter = DecTriMeshVertex<Value>::vv_begin(); iter != DecTriMeshVertex::vv_end(); ++iter)
		{
			//Fix me later. I will define a new class Polygon. 
			center += iter->m_coord ;
			++count ;
		}
		center /= ValueType(count); 
		

		
		return this->m_local_area ; 
	}
	void hello ()
	{
	}
private :
	Value m_local_area ;
};
/**
 * Class of Discrete external triangular mesh face. 
 * Face is assumed to be triangle. 
 */
template <typename Value>
class DecTriMeshFace : public Face<Value>
{
protected :
	typedef typename DecTriMeshPrimitive<Value>::ValueType    ValueType; 
	typedef typename DecTriMeshPrimitive<Value>::FaceType     FaceType;
	typedef typename DecTriMeshPrimitive<Value>::VertexType   VertexType;
	typedef typename DecTriMeshPrimitive<Value>::HalfEdgeType HalfEdgeType ; 
	typedef typename DecTriMeshPrimitive<Value>::EdgeType     EdgeType ; 
	typedef typename DecTriMeshPrimitive<Value>::VectorType   VectorType ;
	typedef typename DecTriMeshPrimitive<Value>::PointType    PointType ;

	typedef typename Face<Value>::FVIter FVIter ;
	typedef DecTriMeshFace<Value> SelfType;


public:

	/**
	 * Calculate the mass center. 
	 */
	void update_mass_center () 
	{
		PointType center (0, 0, 0);
		for (FVIter iter = SelfType::fv_begin(); SelfType::fv_end(); ++iter)
		{
			VertexType * v = (*iter) ;
			center += v->m_coord ;
		}
		center /= 3.0 ;
		this->m_mass_center = center ;
	}

	/**
	 * Calculate the center of circum center. 
	 */
	void update_circum_center () 
	{
		std::vector<PointType *> vset ;
		
		Face2VertexCoord (vset); 
		PointType A =  * vset [0];
		PointType B =  * vset [1];
		PointType C =  * vset [2];

		VectorType pa = A - C; 
		VectorType pb = B - C ;

		VectorType normal = pa ^ pb ;

		VectorType pd1 = pa ^ normal ;
		VectorType pd2 = pb ^ normal ;


		PointType p1 = 0.5 * (A + C);
		PointType p2 = 0.5 * (B + C);

		Line<ValueType> l1 (p1, pd1) ;
		Line<ValueType> l2 (p2, pd2) ;

		this->m_circum_center = l1.intersect(l2); 


//		std::cout<<"-------------"<<m_circum_center<<std::endl;
	}

	/**
	 * 
	 */
	void update_hybrid_center () 
	{
	}



    /**
	 * Get the center of mass ;
	 */
	const PointType & get_mass_center () const 
	{
		return this->m_mass_center ;
	}

    /**
	 * Get the center of the circum-circle; 
	 */
	const PointType & get_circum_center () const
	{
		return this->m_circum_center ;
	}

    /**
	 * Get the hybrid center. 
	 * If the center of the circum-circle is outside the face, we move it to be the center of the edge which opposites the obtuse angle. Otherwise, it is m_circum_center ;
	 */
	const PointType & get_hybrid_center () const 
	{
		
	}
protected :
	
	PointType m_mass_center ;
	PointType m_circum_center ;
	PointType m_hybrid_center ;
};
template <typename PrimitivesType>
class DecTriMesh : public Mesh <PrimitivesType>
{
public :
	DecTriMesh () {}
	virtual ~ DecTriMesh() {}

private :

};


#ifdef SHADE 
class Plane {
public:

    /**
     * Constructs an uninitialised plane.
     */
    INTERFACE Plane() {}

    /**
     * Constructs an initialised plane.
     * @param a The a factor.
     * @param b The b factor.
     * @param c The c factor.
     * @param d The d factor.
     */
    INTERFACE Plane(double _a, double _b, double _c, double _d) : a(_a), b(_b), c(_c), d(_d) {}

    /**
     * Constructs a plane from three points.
     * @param p0 1th point.
     * @param p1 2nd point.
     * @param p2 3th point.
     */
    INTERFACE void fromPoints(const Vector3 &p0, const Vector3 &p1, const Vector3 &p2)
	{
		Vector3 v0 (p0 - p1);
		Vector3 v1 (p2 - p1);
		Vector3 n = v1 ^ v0 ;
		n.normalize();
		a = n.x ;
		b = n.y ;
		c = n.z ;
		d = -(p0.x * a + p0.y * b + p0.z * c);
	}

    /**
     * Constructs a plane from three points and don't normalize the equation.
     * @param p0 1th point.
     * @param p1 2nd point.
     * @param p2 3th point.
     */
    INTERFACE void fromPointsNN(const Vector3 &p0, const Vector3 &p1, const Vector3 &p2)
	{
		Vector3 v0 (p0 - p1);
		Vector3 v1 (p2 - p1);
		Vector3 n = v1 ^ v0 ;
		a = n.x ;
		b = n.y ;
		c = n.z ;
		d = -(p0.x * a + p0.y * b + p0.z * c);          
	}

    /**
     *  Constructs a plane from a point and a normal.
     * @param p The point that lies on the plane.
     * @param n The normal of the plane.
     */
    INTERFACE void fromPointAndNormal(const Vector3 &p0, const Vector3 &n)
	{
		Vector3 normal = n ;
		normal.normalize();
		a = normal.x ;
		b = normal.y ;
		c = normal.z ;
		d = -(p0.x * a + p0.y * b + p0.z * c); 
	}

    /**
     * Makes a dot product between a point and the plane. Result becomes closer to d's factor as the point is near to the plane.
     * @param p The point.
     * @return The dot product.
     */
    double dot(const Vector3 &p) const
    {
		return a * p.x + b * p.y + c * p.z;
    }

    /**
     * Calculates the distance between a point and the plane.
     * @param p The point.
     * @return The distance.
     */
    double dist(const Vector3 &p) const
    {
		return a * p.x + b * p.y + c * p.z + d;
    }

    /**
     * Reflects a vector on the plane.
     * @param p The vector to be reflected.
     * @return The reflected vector.
     */
    Vector3 reflect(const Vector3 &vec)
    {
        double d = dist(vec);
        return  vec + Vector3(-a, -b, -c) * d * 2.0f;
    }

    /**
     * Projects a point onto the plane.
     * @param p The point to be projected.
     * @returned the projected point.
     */
	Vector3 project(const Vector3 &p)
    {
        double h = dist(p);
        return Vector3(   p.x - a * h,
						  p.y - b * h,
						  p.z - c * h);
    }

    /**
     * Tells if a point is on the plane or not.
     * @param p The point to test.
     * @param (optional) The threshold to use to determine if the point is on the plane or not. Indeed we couldn't just test if the distance from the point to the plane is zero because of potential numerical error. Default is 10^-3.
     * @return true if the point is on the plane, false otherwhise.
     */
    bool isOnPlane(const Vector3 &p, double threshold = 0.001f)
    {
        double d = dist(p);
        if (d < threshold && d > -threshold)
            return true;
        return false;
    }

    /**
     * Determines if a line and the place intersect.
     * @param p0 The first point of the line.
     * @param p1 The second point of the line.
     * @param t A value that will be filled with the distance to travel to get the intersection.
     * @return If the plane and the line are parallel false is returned, true otherwhise since the line obviously intersects with the plane.
     */
    bool intersectWithLine(const Vector3 &p0, const Vector3 &p1, float &t)
    {
        Vector3 dir = p1 - p0;
        double div = dot(dir);
        if (div == 0)
            return false;

		t = -dist(p0) / div;
        return true;
    }

private :
    /**The a factor of the plane.*/
    float a;
    /**The b factor of the plane.*/
    float b;
    /**The c factor of the plane.*/
    float c;
    /**The d factor of the plane.*/
    float d;
};


/**
 * Determine the gradients of function f on mesh. 
 * @param P[3] are the three points in one triangle face,
 * @param f[3] are the corresponding value on 3 vertice of function f. 
 */
void gradient (Point p[3], float f[3])
{
	Plane plane; 
	plane.fromPoints(p[0], p[1], p[2]); 
	// Construct the first base function .
	Vector3 v1 = p[0] - p[2]; 
	
	// Construct the second base function .
	
}

#endif

#endif
